Wilson Fermion Determinant in Lattice QCD 



Keitaro Nagata 
Department of Physics, The University of Tokyo, 
Bunkyo-ku, Tokyo 113-0033 JAPAN 

Atsushi Nakamura 

Research Institute for Information Science and Education, 
Hiroshima University, 
Higashi-Hiroshima 739-8527 JAPAN 
(Dated: September 14, 2010) 

We present a formula for reducing the rank of Wilson fermions from 4N c N x N y N z Nt to 4N c N x N y N z keep- 
ing the value of its determinant. We analyse eigenvalues of a reduced matrix and coefficients C„ in the fugacity 
expansion of the fermion determinant 2~2 n C n (exp(/i/T)) n , which play an important role in the canonical for- 
mulation, using lattice QCD configurations on a 4 4 lattice. Numerically, log \C n \ varies as N x N y N z , and goes 
easily over the standard numerical range; We give a simple cure for that. The phase of C'„ correlates with the 
distribution of the Polyakov loop in the complex plain. These results lay the groundwork for future finite density 
calculations in lattice QCD. 
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L INTRODUCTION 

QCD at finite temperature and density has been one of the most attracting subjects in physics. Many phenomeno- 
logical models predict that the QCD phase diagram is expected to have a very rich structure, and thoroughgoing 
analyses of heavy ion data have been made to show that we are sweeping finite temperature and density regions. 
See Ref. fl]. 

First-principle calculations based on QCD are now highly called. If such calculations would be at our hand, 
their outcomes are also very valuable for many research fields: high energy heavy ion collisions, the high density 
interior of neutron stars and the last stages of the star evolution. Needless to say, the inside of nucleus is also a 
baryon rich environment, and lots of contributions to nuclear physics could be expected. 

Unfortunately, the first principle lattice QCD simulation suffers from the sign problem. Nevertheless, there have 
been many progresses such as the reweighting method fl, the imaginary chemical potentially, Q] and the canonical 
formulation g|; now some light is shed on the QCD phase diagram. 

Most of lattice QCD studies with non-zero density were done with the use of staggered fermions. It is desirable 
to study lattice QCD with Wilson fermions because it is free from the fourth-root problem. At zero density, thanks 
to several algorithm developments, lattice QCD simulations with Wilson fermions are now possible even on the 
physical quark masses. 

In case of the lattice QCD simulations with finite chemical potential fi, often we must handle the fermion 
determinant dct A(/i), directly. For example, the reweighting method requires a ratio of two determinants, 

detAQ/) 

det AO) ■ 1 ; 



The density of state method needs the phase information fl. The canonical formulation needs the Fourier transfor- 
mation of the fermion determinant 

det A n = i- J d (^) e~ m ^l T det A( W ), (2) 

with the quark number n and the imaginary chemical potential /i/. In these approaches, the heaviest part of the 
numerical calculations is the evaluation of the determinant. An efficient way of the determinant evaluation is highly 
desirable. It is very useful if we can transform the fermion matrix A into a compressed one whose rank is less 
than the original one, and yet it gives the same value of the determinant, since the numerical cost to evaluate a 
determinant is usually proportional to the third power of the matrix rank. 

Such a transformation was found for the staggered fermion by Gibbs |8] and Hasenfratz and Toussaint fl, and 



used in finite density simulations, e.g. II 1 CM — I- 1 3fl . Their method has also an advantage in the canonical formulation. 



With the reduction method, the fermion determinant is expressed in powers of fugacity, 

dctA( A1 ) = ^a i (e' l/T ) . (3) 

n 

If we obtain the coefficients C„ , the Fourier transformation in the canonical formulation is easily carried out. 

A reduction method for Wilson fermions has not been established yet. It is unfeasible to apply the method 
for staggered fermions in Isl to Wilson fermions in a naive way because of singular parts contained in the 
Wilson fermion matrix. Expansions based on the trace-log formula have been proposed for the Wilson fermion 



determinant 1 141 - 121 II . An efficient method to calculate exactly the Wilson fermion determinant is valuable for finite 



density simulations with Wilson fermions. 

The purpose of the present work is to construct a reduction method for Wilson fermions. In Ref. Il22ll . Borici 
derived a reduction method that can be applied to Wilson fermions, and tested it using a Schwinger model (QED2) 
with staggered fermions. We develop further the method of \l2^ and derive a reduction formula, which rearranges 
the Wilson fermion determinant in powers of fugacity and reduces the numerical cost. 

Similar to the method in J^Hgl, the Wilson fermion matrix is expressed in a time-plane block matrix form. 
Projection operators contained in the Wilson fermion matrix make it possible to transform forward and backward 
hopping parts separately. Owing to the property of the projection operators, the Wilson fermion matrix is trans- 
formed so that the determinant in the time-plane block form can be carried out analytically. The determinant of 
the Wilson fermion is then reduced into that of a reduced matrix, whose size is smaller than the original one. The 
problem results in the diagonalization of the reduced matrix instead of the original matrix. Solving the eigenvalue 
problem for the reduced matrix, the Wilson fermion determinant is expressed in powers of fugacity. 

This paper is organized as follow. In the next section, we show the reduction method for the Wilson fermions. 
In section [III] as an illustration, we perform numerical simulations on a small 4 4 lattice and calculate the Wilson 
fermion determinant using the reduction method. We discuss the properties of the coefficients of the fugacity 
expansion. The results are not to be regarded as physical, due to the small lattice size, but lay the groundwork for 
future realistic calculations. The final section is devoted to a summary. In the appendix, we give (1) the detail of 
the calculation of the determinant of a permutation matrix P used in the reduction formula, (2) a simple numerical 
trick to evaluate the fugacity expansion coefficients, and (3) a possible alternative formulation. 
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II. FRAMEWORK 



A. Structure of Fermion Matrix 



We employ the Wilson fermions defined by 

3 



A(x,x') = 5 x , x , - K^{(r-7i)^(^ x , ;X+ i + (r + 7i )^(^x',x-i} 



i=l 



+ Sci over ? 

S Clover = — Sx.x'CswK ^ V(j,yF(xv' 



(4) 



where r, k and [i are the Wilson term, hopping parameter and chemical potential, respectively. We include the 
clover term with the coefficient Csw- For later convenience, we divide the quark matrix into three terms according 
to their time dependence 



A = B - 2z- 1 nr-V - 2znr+V f . 
Here r± = (r ± j4)/2 and z = e~ M , and 

3 

B(x,x') ee 6 X , X , -«^{(r- li)Ui(x)5 x ,^ +l + (r + li)Uj (x')6 x ,^} 
i=i 

+ Sci over •> 

V(x,x') ee Ui{x)5 x , x+h 



(5) 



(6) 
(7) 
(8) 



They satisfy = I. Note that r± are projection operators in the case that r = 1. In a time-plane block matrix 
form, B and V are given by 



t'=l 



t'=N t 



t = 1 
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B. Reduction Formula for Wilson Fermions 



Now, we derive a reduction formula for the Wilson fermions. A starting point is to define a matrix B22I1 . 

P={c a r_+c b r+Vz- 1 ), (10) 

which is referred to as a permutation matrix [ 220 . The parameters c a and c b are arbitrary scalar except for zero, 
and may be set to one. We can use these parameters to check the following reduction formula numerically. Since 
r± are singular, the matrix P must contain both of them; otherwise P is singular. It is straightforward to check 
det(P) = (c a c b z~ 1 ) N/2 , where N = 4N c N x N y N z N t . Multiplied by P, the quark matrix is transformed into 



AP = (c a Br- - 2c b nr + ) + (c b Br+ - 2c a nr-)Vz^ 1 . 
In the time-plane block matrix form, the first and second terms of Eq. (fTTT i are given by 

(c a Br_ - 2c b nr + ) = 



a 2 



V 



a Nt j 



( fcz- 1 



(c b Br + — 2c a Kr J)V z 



foz- 1 




/ 



The block-matrices are given by 



a i = a ab ^(x,y,t i ) 

= c a B ab ^{x, y, U) r™ - 2c b n rf 6 ab 6(x - y), 

= c b B ac ^(x, y, U) rl v Uf{y, t t ) - 2c a K r^6(x ~ y)Uf(y, t t ), 



(11) 



(12) 



(13) 



(14) 
(15) 
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where the dimensions of a,; and /3{ are given by iV re< i = N/Nt = 4N x N y N z N c . We factor out a negative sign 
caused by anti-periodic boundary conditions from the definition of f3pf t . Therefore the negative sign appears at the 
lower-left corner in Eq. dl3) . The two block-matrices have different meaning; contains only spatial hopping 
terms with a fixed time t = U, while /3, contains temporal hopping terms as well as spatial ones due to temporal 
link variables. 

Combining Eqs. ( fT2b and JT3b , we can carry out the determinant in the time-plane block matrix 



det AP = 



Oil 



a 2 faz~ x 



\ 



o 3 



JJdet(tti) det (1 + z- Nt Q) 

\i=l J 



(16) 



where Q — {a x 1 j3\) ■ ■ ■ [a N j9jv«)> which we refer to as a reduced matrix. Substituting det(P) = (c a ci,z l ) N / 2 , 
we obtain 



det A = {c a c b )- N / 2 z- N/2 ( JJdet(oi) J det (z Nt + Q) 



(17) 



Here, the rank of the matrices «j and Q is given by iV rc d = N/Nt, while that of the Wilson fermion is originally 
given by N. Thus the reduction formula makes the computation of the determinant 1 /Nf less time. Furthermore, 
the /i dependent parts are separated from the hopping terms, and appear at the overall factor and the second 
determinant. 
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FIG. 1: The left panel depicts closed loops on a plane t = t%, which contribute to det on, while 
the right one does paths of quarks from t = t% to t = ijv t , which contribute to the reduced 
matrix Q. 



Equation ( fTTI i consists of sub-determinants: J^det(ai) and det(z Nt + Q). As we have explained, oij describes 
spatial hopping terms at a time-slice with 4j, Hence, det ai describes closed loops in a plane with a fixed time 
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FIG. 2: Schematic figures for the reduction procedure. 

t = U (left panel in Fig.Q}. On the other hand, Pi contains temporal hopping terms from t = t j to t = The 
reduced matrix Q describes paths of the propagation of quarks from t = t\ to t = t^ t (right panel in Fig. [TJ. Thus, 
the reduction formula separates closed loops at time-slices from paths of the propagation of quarks (Fig. [2j. 

Now, we solve an eigenvalue problem det(Q — XI) = 0. With the eigenvalues A, the determinant of the reduced 
matrix is written as 

AW 

det(z N * +Q) = Y[(\ n + z N <), (18) 

71=1 

which is expanded in powers of z Nt , 

JV rod AW/2 
n=l n=-N [cd /2 

= C - Wrod/2 (^ t ) JV " d/2 + • • • + c + • • • + c Nnd/2 (z N T N - d/2 , d9) 

where we replace c n by c_„ to obtain the second line from the first one. This is an expansion with regard to 
(inverse) fugacity z * — exp(— n/T). Equivalently, this can be interpreted as a winding number expansion, 
because z Nt comes from closed loops that make a round the lattice in the time-direction. Note that the expansion 
Eq. (fT9l is exactly done and does not involve any approximation. 
Finally, we obtain the reduced quark determinant 

AW 2 

detA(/i)= C ^ /T ) n , (20) 

n =-A r red/2 

Here, C n = Cc n with C = (c a c b y N / 2 (n^i det(oi)). 

The coefficients c„ have two properties. If a chemical potential is pure imaginary \i = i^j, then (z)* = z^ 1 and 
(det A(//r))* = (det(A (///))). These conditions bring about the first property c* = c_ n . Note that C-N lad /2 — 

CjV red /2 = ]ln=l = 1- 

The second property is concerned with the center transformation Z3. Under Z3 transformation, the time com- 
ponents of the link variables are transformed as 



U 4 (ti) wU 4 (U), 



(21) 
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where w — exp(2iri/3) is an element of Z3. Regarding the n-th term in Eq. (120b . if the winding number n is a 
multiple of N c , the coefficient c„ is Z3 invariant, otherwise c„ is not Z3 invariant. Thus, c n are classified in terms 
ofZ 3 

{center invariant (n = 3m) 
(22) 
center variant (71 = 3m + 1, 3m + 2) 

where m is an integer. It is known that the center symmetry is explicitly broken in the presence of quarks. In the 
quark determinant, the explicit breaking of the center symmetry is caused by the terms having winding numbers 
not multiple of N c . 

III. NUMERICAL RESULTS 

In this section, we demonstrate the calculation of the quark determinant dot A(/i) using the reduction for- 
mula. In order to see the temperature dependence of detA(/x), we set /3 — 1.85 and 2.0. We employ 
(k,C sw ) = (0.14007,1.5759) and (0.1369,1.5058) for (3 = 1.85 and 2.0, respectively. We perform hybrid 
Monte Carlo (HMC) simulations on the 4 4 lattice with 1,000 quench updates and 100 full QCD HMC trajectories 
as thermalization. After the thermalization, we measure the quark determinant on four configurations separated 
by 20 HMC trajectories between measurements. Fundamental numerical data of the configurations used for the 
measurements are shown in Table HI and HI1 

TABLE I: The values of the determinant with fj, — , Polyakov loop and plaquette for /3 = 1.85. 
(i), (ii), (iii) and (iv) correspond to the configurations measured. 

det A(0) Polyakov loop Plaquette 

(i) 3.0957 x 1(T 19 0.04377 - 0.25418i 0.53338 

(ii) 2.0921 x 10~ 21 -0.03234 + 0.0871H 0.50668 

(iii) 2.2560 x 10~ 21 -0.16365 - 0.10135i 0.52471 

(iv) 5.1115 x 10 -18 0.49234 - 0.12163i 0.53313 



TABLE II: The values of the determinant with /1 = 0, Polyakov loop and plaquette for ft — 2.0. 



detA(0) Polyakov loop Plaquette 

(i) 8.7586 x 10~ 12 0.37590 + 0.004H 0.57810 

(ii) 3.0329 x 10~ 12 0.13827 - 0.1978i 0.57107 

(iii) 1.1159 x 10" 12 -0.22324 - 0.4285i 0.57491 

(iv) 1.2578 x 10~ 12 -0.35711 - 0.6028i 0.57954 



One of the advantages of the reduction method is that it makes easy to calculate the fi dependence of the quark 
determinant. Once we perform the reduction procedure and obtain A„ or c„, we can obtain det A(/i) for arbitrary 
/i. In Figs.[3]and[4] we show the /1 dependence of the determinant. The values remain near the starting points when 
fj, is small, and move rapidly when fj, exceeds around 0.9. 



8 





5e 


014 











-5e 


014 


< 






o5 


-1e 


013 


Q 






E 


-1.5e 


013 




-2e 


013 




-2.5e 


013 



-5e-014 





5e-022 









-5e-022 


< 




CD 


-1e-021 


Q 




E 


-1.5e-021 




-2e-021 




-2.5e-021 




5e-014 
Re[Det A(n)] 



1e-013 




-4e-021-2e-021 2e-021 4e-021 
Re[Det A(n)] 



CD 

Q 



CD 

Q 



5e-019 


-5e-019 
-1e-018 
-1.5e-018 
-2e-018 



-2.5e-018 



-2e-018 



140 
120 
100 
80 
60 
40 
20 

-20 



(ii) 





Re[Det A(n)] 



2e-018 



-100 



(iv) 



100 200 300 400 500 
Re[Det A(n)] 



FIG. 3: Parametric plot of Det A(/i) from \i = to [i = 1 for /3 = 1.85. The points are denoted 
for Sfi = 0.01. 



Next, we study the reduction method in more detail. The distribution of the eigenvalues A of the reduced 
matrix Q is shown in Figs. [5] We observe that the eigenvalues are split in two regions. Almost half of the 
eigenvalues are distributed in a region |A| > 5, and the other half in a region |A| < 0.5. There is a margin 
between two regions, where no eigenvalue is found, as we can see in the right panels in Fig. [5] The splitting of 
the eigenvalues is observed in the eight measurements. Note that the eigenvalues are constrained by the condition 
11^=1 = 1- Qualitatively, this can be understood from the fact that the matrix Q is a product of the block 
matrices A; = (a" 1 /^). It is expected that when the system is in equilibrium A; moderately depends on time. In 
such a case, we can express A ~ A + SAi, where A is independent of time. Assuming the time-dependent part 
8Ai is small, Q = Ui=i A % ~ A Nt + 0(6). Then, A Nt causes the splitting of the eigenvalues of the matrix Q for 
eigen(A) > 1 and eigen(A) < 1 cases. 

The coefficient c n is a polynomial of the eigenvalues A according to Eq. ( [T9l l. Because the number of the 
eigenvalues iV re d is large, there appear two numerical problems. First problem is for an accuracy. We employ a 
recursive method in order to determine c„ in enough precision. Second problem is that c„ exceeds the range where 
a number can be represented in double precision: about 10~ 308 ~ 10 308 . In order to overcome this problem, we 
develop a special routine to extend exponential part. See Appendix|B] For the check, we compare our results with 
those obtained by using FM multi-precision library (FMLIB) ll23ll . 

We plot the absolute value of c„ as a function of the winding number n in Figs. [6] where we show the results 
only for the configurations (i) both in high and low temperatures because results for the other configurations are 
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FIG. 4: Parametric plot of Det A(/i) from /j, = to fj. = 1 for ft = 2.0. The points are denoted 
for 5 ft = 0.01. 



very similar to Figs. [6] As we have mentioned, |c„| goes over the standard numerical range and reaches at 10 900 
at most, which is much larger than the maximum value in double precision. Note that the overall factor C is order 
10 -900 , then the cancellation between C and c n makes their product C n — Cc n ordinary order. For both /3 = 1.8 
and 2.0, we find that \c n \ is maximum at n = and decreases exponentially as \n\ becomes larger. 

Next, we show the absolute value of C n e ntl / T for several chemical potentials in Figs. UJ and [8] In contrast to 
|c n |, the fugacity factor e n ^l T becomes larger as n becomes larger. The difference between the n-dependence of 
|c n | and e n ^l T leads to a peak for \C n e ntJ, ^ T \, as we can see in Figs. [7] and [8] Several terms in the vicinity of the 
peak dominate det A(/i). For instance, det A(0) is dominated by terms near n = 0. The location of the peak 
moves towards larger values of n as fi becomes larger. However the /i dependence of the location of the peak is not 
so strong. Even for the chemical potential near to fi = 1, significant contributions come from terms with n < 100. 
In the following, we consider terms with n < 100. 

The phase of c„ as a function of n are shown in Figs.|9]and Figs.QI)] In all the eight configurations, there is a 
symmetry under a rotation with 7r. This symmetry indicates the relation c* = c_ n is numerically satisfied. The 
phase of c n complicatedly depends on the winding number, temperature and configuration. We find that there are 
two particular n-dependence. One is that the phase of c„ is a continuous function of n, e.g. see (iv) in Figs. [9] (This 
means that we can fit the phase of c n in terms of a continuous function, although n is not a continuous variable.) 
The other is that the phase of c n is spit into three lines classified in terms of mod(n, 3), e.g. see (i) in Figs. [9] Each 
line is a continuous function of n = 3m, 3m + 1 or 3m + 2 and there is a gap about 27r/3 between lines. This 
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FIG. 5: The distribution of the eigenvalues A in the complex plane. The topped and bottomed 
panels are for /3 = 1.85 and 2.0, respectively. The left and right panels show the distributions 
in two different scales. 




splitting causes the cancellation of C n e nii l T among neighboring three terms with n = 3m, 3m + 1, 3m + 2 and 
suppresses the magnitude of the determinant. For instance, the four values of det A(/i) corresponding to the four 
configurations in f3 = 1.85 are classified in two groups, as we have seen in Figs.[3j det A(/i) is of order 10~ 20 in 
the configurations (i), (ii) and (iii), and it is of order 100 in the configuration (iv). This observation is related to the 
behavior of the phase of c n . 
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FIG. 7: The distribution of |C„e nM/ ' T | for /3 = 1.85. The four panels correspond to the config- 
urations (i), (ii), (iii) and (iv), respectively. 



IV. SUMMARY 

In this paper, we have presented the reduction formula for the Wilson fermion determinant. The formula reduces 
the numerical cost to evaluate the Wilson fermion determinant. The point is that the Wilson fermion matrix contains 
the projection operators, which enable to transform the fermion matrix so that the temporal part of the determinant 
can be performed analytically. Thus the Wilson fermion determinant is reduced to the determinant of the reduced 
matrix. Solving the eigenvalue problem for the reduced matrix, the determinant is expressed in powers of fugacity. 
Although the basic idea for the reduction method is similar to that for staggered fermions, a difference comes from 
the use of the projection operators. 

We perform the numerical simulations on the 4 4 lattice and calculate the Wilson fermion determinant using the 
reduction formula. In order to determine the coefficients of the fugacity expansion in enough accuracy, we employ 
the recursive method and develop the special routine. Furthermore, we compare our results with those obtained by 
using a multi-precision library. 

We discussed the properties of the eigenvalues of the reduced matrix and of the coefficients c„ of the fugacity 
expansion. The eigenvalues show an interesting behaviour; they are split in two regions. We find that there are two 
particular behaviours for the winding number dependence of the phase of the coefficients. One is that the phase 
of c„ is a continuous function of the winding number. The other is that the phase of c„ is split into three lines 
classified in terms of mod(n, 3) with the gap about 2tt/3 between lines. 
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FIG. 8: The distribution of |C"„e nfl/T | for f3 = 2.0. The four 
panels correspond to the configurations (i), (ii), (iii) and (iv), 
respectively. 
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Appendix A: Determinant of permutation matrix 



Here we calculate the determinant of the permutation matrix P. Since the projection operators are singular 
det(r± ) = 0, we need to use a simple trick in order to obtain det P; first we reduce the determinant in the case that 
r 7^ 1. Then we take the limit r — Y 1, after eliminating the singularity. To perform this, we summarize identities 
of the projection operators for arbitrary r. They are defined by 

r ±74 



Using the definitions, it is straightforward to obtain 



r 2 - 1 



These equations lead to the inverse matrices 



Using Eqs. iA2i . we obtain 



detP = 



(r+) 1 = -r_ , (r_) 1 = -r_, 
e e 



c b r + z 1 Vi ) 2 

c a r_ c h r + z _1 V2,3 
c„r_ 



(Al) 



(A2a) 



(A2b) 



= det fcfr^+cfV^^n^A+i 
Considering the Dirac components, (r±) Nt are given by 



Cbr+z 1 W«-l,iV t 



(A3) 



/m 2 



(A4a) 



(r_) 



v 



\ 



m 2 / 



(A4b) 



Having these two terms, there is no singularity in Eq. JA3b . Then, we obtain 

detP = z- N / 2 (c a c b ) N / 2 . 



(A5) 



Appendix B: Calculation of coefficients c„ 

As we have shown, the coefficients c„ in Eq. ( fT9b vary from order one to order 10 900 even on the small 4 4 lattice. 
They cannot be handled in the double precision. This problem usually happens when we consider expansions of the 
fermion determinant. So far, arbitrary accuracy libraries are often employed in order to calculate the coefficients. 
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We calculate c n as follows in a recursive way: 

M M-1 



c k e = (B Q + srf) J2 c ^ k (Bi) 



k=0 k=0 



and 



Cq — BqCq (B2) 
C£ = Bk-iCk + B k C k -! (k = l,2,--- ,M-1) (B3) 
= BiCm-i (B4) 

In order to express Cfc, we need wide range of floating numbers, but in Eqs. dB3b we do not need very high 
precision. In other words, we need wide range of the exponent, but we do not need very large significant numbers. 
We express each real and imaginary parts of C k in a form of 

a x b L (B5) 

where 

1 < \a\ < b (B6) 

and a is a double precision real and L is an integer. When we solve the recursion relation Eqs. ( IB3b . we express all 
Cfc, C' k and B k in this form. The base b can be any number, and we set it to be 8. 

To see if this simple trick works or not, we calculate several cases by this method, and by a high accuracy 
Ubrarv, FMLIbB. We got the same results. Although this method works for obtaining the coefficients, Cj in a 
sufficient double precision, we found a peculiar configuration on which a huge cancellation occurs in the sum of 
Cj x exp(ifj,/T) and the double precision is not enough to get a correct value of the determinant. 



Appendix C: Alternative approach 

In this appendix, we give another possible transformation of the Wilson fermion determinant. It is more direct 
extension of the Gibbs's approach for the staggered fermion, and may give a general base. Unfortunately, in 
present-days numerical algorithms we cannot find a reliable one to solve a generalized eigenvalue problem if 
involved matrices are singular. But if in future this problem is solved, the following can be another good starting 
point. 

Keeping in mind that for Wilson fermions, {—k(t — J^V)^ 1 does not exist unless the Wilson term r ^ 1, we 
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can apply similar transformation in the staggered fermion case by Gibbs to the Wilson fermion, and get 



detA = dcti(zB + z 2 (-K(r + 7 4 )y t ) + (-K(r-7 4 )K)) 

= dct - z (zBV + z 2 (- K (r + 74)) + (-/s(r - 1a)V 2 )) V' 1 



= z 



-N 



I dot V 



= z- N 

Here the block-matrices are given by 

BV = 

( 



BV — z(—k(t + 74)) I 
n(r — "n)V 2 —z 

I -n(r + 74) 



-BV I 

n(r-j 4 )V 2 










BiV 12 














B2V23 ■ ■ ■ 














• BjV t -2Viv t -2JVt-l 














BN t -\V Nt -i Nt 














and 



V 









V12V23 



















V23V34 • • • 





















VN t -2 N t -1 VN t -INt 


VjV t -lJV t VAT t l 


















Vjv t lVi2 














(CI) 
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By exchange columns and raws, 

-BV I 

n{r- l4 )V 2 







-SiVi 2 1 
aV Nll Vi2 


















—B2V23 1 
aVi 2 V 2 3 
























— BN t -lVN t -lN t 1 

aVlv t -2Jv t -iVjv t -iJv t 


— B N t VN t i 1 
o:Viv t -iJv t VAr t i 















Here we introduce a = n(r — 74). Applying the same exchange of the columns and raws, 



-n(r + 74) 
I 



-n(r + 74) 
1 


















— K(r + 7 4 ) 
1 




























-«(r + 74) 
1 



(C2) 



(C3) 



Then we can write 



where 



T = 



det A 


= z- N det(T- 


zS) 


( ° 


h 


••• 


\ 








t 2 ■■■ 











••• 








• • • t Nt - 2 














tN t -l 










/ 


k = ^ 




-BiV iti+ i 1 \ 






- 74)^-1,1^+1 / 



(C4) 



(C5) 
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and 



s 





••• 








s 


••■ 











••■ 


















s 











... o 


s 




(C6) 



Each U and s is (4N c N x N y N z ) x (4N c N x N y N z ) matrix. 

Eq. ( IC4I ) is a form of the generalized eigenvalue problem[24]. There is a mathematical theorem (Generalized 
Schur Decomposition) which tells us that there exist unitary matrices Y and Z such that Y' SZ and Y^TZ are 
upper triangular. Let a k and /3k be diagonal elements of these matrices. Then 



det(T - zS) = det FZ* JJ(a fc - zfi k ) 

k 



(C7) 



Half of a's and /3's vanish. 

This formula has an advantage that T and S do not have inverse matrix like Q in Eq. (fTTI ). and can be easily 
constructed. A problem is that matrices T and S are singular. To our knowledge, no stable algorithm is know to 
solve the generalized eigenvalue problem in such a case. 
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